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In this paper, we use the Quantum Trajectory Method, a Monte Carlo Simulation applied to 
quantum systems, to study a microlaser operating on three-level atoms interacting with a two-mode 
cavity. We are interested in the quantum statistical properties of the cavity field at steady state. 
It's shown that steady state solution does exist when the detailed balance condition doesn't apply. 
We highlight the diff'erences between a single mode microlaser and a two-mode microlaser. Special 
attention is paid to the one-photon trapping state. The second order correlation function g^^\r) of 
a single mode is studied. We show the effects of the velocity spread of the atoms used to pump the 
microlaser cavity on the second order correlation function, trapping states and phase transition of 
the cavity field. We illustrate an interesting effect of the atomic velocity spread on the coherence 
function when the cavity field of a mono- velocity atomic beam microlaser exhibits anti-bunching. 



I. INTRODUCTION 

The single atom laser (microlaser) has been used suc- 
cessfully in the last decade as a source of non-classical 
light [H m [3] . Following the same principles of the sin- 
gle atom maser (micromaser) , the microlaser operates by 
pumping a high finesse optical cavity by a beam of ex- 
cited atoms. The microlaser we are considering in this 
paper is different than the single atom laser first operated 
by McKeever where one and the same atom is used. 
We use the "microlaser" term in this paper however to 
describe the experimental setup of the MIT microlaser 
[T]. Being one of the best systems to study cavity QED 
effects, the single atom maser /laser has received a lot of 
interest in the quantum optics community. One interest- 
ing type of microlaser utilizes bi-modal cavities pumped 
by multi-level atoms |5J [SJ [7]. A two- mode micromaser 
pumped by three-level atoms in A-configuration has been 
analyzed analytically by F. L. Kien et. al. in [7 . In the 
current work we use the Quantum Trajectory Method 
[HJinilinjIII] to analyze the two-mode microlaser operated 
by pumping a doubly resonant optical cavity by atoms 
characterized by two lasing transitions and strongly cou- 
pled to the cavity. Atoms used to pump the microlaser 
cavity are produced by an oven and the velocities of the 
atoms obey a thermal velocity distribution. The atoms 
then are let to pass through a velocity selector to unify 
the speeds of the atoms and consequently their interac- 
tion times with the cavity. Since the efficiency of the 
velocity selector is not perfect, atoms passing through 
the cavity still have a slight velocity variation. Our main 
target in this paper is to calculate the cavity field coher- 
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ence function and investigate the effects of the variation 
of the atoms' velocities on the coherence function and 
other statistical properties of the cavity field. This pa- 
per is divided into five sections. In Sec. II we review 
the theory of the two-mode microlaser and the coherence 
function of the field inside the microlaser cavity and elab- 
orate on the interesting one-photon trapping state. In 
Sec. Ill we introduce the Quantum Trajectory Method 
and apply it to the two-mode microlaser and show the 
results of calculating the coherence function ^'^•'(t) nu- 
merically and analytically. In Sec. IV, we show the re- 
sults of including the variation of the atoms' velocities in 
the numerical simulation and elucidate its effect on the 
statistical properties of the microlaser field. We conclude 
our work in Sec. V and present suggestions for further 
investigations. 

II. THEORY 

A schematic diagram of the energy levels of the three- 
level atoms used to pump the bi-modal microlaser cavity 
is shown in Fig. [T] In the most general case, each mode 
has its own angular frequency tOa , coupling strength with 
the atom ga, decay rate ja, and mean number of ther- 
mal photons rif,^ where a = 1,2. Although, unlike the 
micromaser, the mean number of thermal photons for a 
microlaser is typically zero, we include nbi,nb2 in this pa- 
per for the sake of generality. All the atoms are excited 
to the higher level \a) before they enter the cavity. Atoms 
are statistically independent (have random arrival times) 
and their dwelling time inside the cavity is much shorter 
than the mean inter-atomic arrival times. 

We assume that the detuning between the two mode 
frequencies is large, compared to the atom- field coupling 
strengths (71,(72, so that each mode interacts only with 
the respective atomic transition. Since the cavity is sub- 
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FIG. 1: A schematic diagram of the energy levels of a three 
level atom in A-configuration 

jected to the incoherent field decay, we are going to use 
density matrix formalism where the field is represented 
by the density matrix p. The diagonal elements of the 



field density matrix P„i.„2, which represent the joint 
probability to find ni photons in mode 1 and n2 photons 
j.n mode 2, constitute a two dimensional matrix. In what 
follows we are going to summarize the procedure given in 
[7] to calculate the steady state solution of Pni,n2 ■ We are 
not interested in the off-diagonal elements of the density 
matrix in the current work. Changes to the density ma- 
trix are induced by two different and independent proce- 
dures, the atom field interaction and the field dissipation 
from the cavity. The change to the density matrix due to 
the interaction between the field and a single atom dur- 
ing the interaction time Tint is represented by the super- 
operator F(t) such as p{ti + Tint) = F (Tint) p{ti) ■ The 
field dissipation is controlled by the LiouviUian super- 
operator Lq given by |12j 



p = Lop 
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where ya[,aij and ya2Ta2j are the field operators of 

the two modes and nbi,nh2 are the mean number of ther- 
mal photons in the two modes. It can be shown [13^ that 
for the random arrival of atoms, which is a Poisson pro- 
cess, the master equation controlling the change of the 
density matrix p is given by 

p^R[F{Tint)-l]p+Lop (2) 

where R is the rate of the atoms injection. At steady 
state, we have 

P = (3) 



exactly compensated by the change due to field dissipa- 
tion from the cavity. In other words the net change per 
unit time of the field density matrix is zero. To obtain 
the form of F (r) , we start by writing the interaction 
Hamiltonian of the atom-field system during the interac- 
tion assuming that the transition between the two ground 
levels is forbidden: 

^^int = T^gi {a\ai + aicrl^ + %2 {alcr2 + a2crj^ (5) 

The Shrodinger's equation in the interaction picture, gov- 
erning the evolution of the atom-field wavefunction is 



R[l-Fin^,)]p^Lop (4) = (6) 

ot 

which simply means that, the change per unit time in the 

field density matrix due to the atom field interaction is where \ip{t)) is generally expressed as: 



IV'(^)) = y2 Ca,n,ni \a, 11, m) + Cb^^n^m |&1, TTl) + Cb^^n.m 1^2,^, w) 



n.m 



(7) 



and n, m are the numbers of photons in mode 1 and and hence 
mode 2 respectively. Since at the beginning of the inter- „ /'^^ _ i |2 (q\ 

action, the atom is in its excited state, then ^nM'J) - |Ca,«,m| {^) 

Inserting the expression of \ip{t)) into the Schrodinger's 

equation we obtain 

l'/'(0)) = ^Ca,„,m |a,n,m) (8) 

,n-\-l.m 

(10) 
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from which we obtain 



(11) 



(12) 



(13) 



By solving ( 11|12|13" ) and taking the origin of time to 
be at t = we get 



Ca,n.m(Tint) = Ca.n.m{0)cOS^J (n + 1) + .g| (to + l)Tint 

(14) 



Cf,i,n,m(Tint) = Ca,n-l,m(0)gi 



Cf,2,„,m(Tint) = CQ,„,,„_i(0)g2V^ 



Defining 



sin^/ gl (n) + g| (to +T)Tint 

(«) + 5l (™ + 1) 
(15) 

sinyfgl (n + 1) + g| (TO)Tint 

VS? ("+l) + .9l(™) 
(16) 



A (n, m) = ^gl {n+\) + gl (to+1) (17) 



we conclude that 



fn,m(Tint) = -P«,m(0)cOS^ [A (n, to) Tint] + 5inP„_i,m(0) 



+g2TOP„,„_l(0) 



sin^ [A (n- l,TO)Tint] 
A^ (n — 1, to) 

sin^ [A {n, to - 1) Tjnt] 
A^ (rt, m — 1) 



We now have the form of F (t) . By plugging it into ([3|, the density matrix, we get: 
Q and rewriting it in terms of the diagonal elements of 



hi N n or ■ 2rw ^ 1 D/ N , 2 sm A(ni - 1, ^nt 

P(ni,n2) = = i?^{-sm [A(ni, r7,2)TintJ P(?^l, ?^2) + r^T ; ^ 

V{ni - 1, 712) 



F(ni - 1,712) 
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sin^ [A(ni,ri2 - l)T{^t] 



P(71l,7l2-1)}- 



A2(?zi,n2 - 1) 

7i("-fci + 1) [("1 + l)^'("-i + 1, "-2) - niP(ni, 712)] + 717161 [7iiP(7ii - 1, 772) - (771 + l)P(77i, 772)] 



(18) 



The previous equation contains the probability flow 
terms between any two successive energy levels cavity 
field as shown in Fig. [2j 

The arrows represent the probability flow between ad- 
jacent energy states. Since at steady state, the net prob- 
ability flow from level (0,0) is zero, we conclude that the 
sum of the probability flow terms A and B is zero. Under 
the symmetric operation of the microlaser, deflned by, 

ffi = .92, 71 = 72, "61 = 



we can conclude that the probability flow bundles A and 
B are equal and hence each of them will be identically 
zero. By induction, one can conclude that the three-term 
probability flow between any two levels in ( 18 ) is zero and 
we obtain |7]: 



2 Sm A(77i-l,772)Ti„t , / , 1\r D^- M \ Ut 1 M 

Rg 771 r^T --, ^ Pyni - 1,772) = l{nb + 1) 77iP(77i, 772) - ini, 77iP(77i - 1,772) 

A^(77l - 1,772) 



(19) 
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FIG. 2: A schematic diagram of the lowest energy states of 
the field inside a bi-modal cavity. The arrows represent the 
probability flow between adjacent energy states. 
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A2(ni,n2 - 1) 
which can be solved together to get 
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P(ni,n2) = P(0,0) 
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where P{0, 0) is obtained from the normahzation 
condition ^ P(7ii,7i2) = 1- The probabihty to find n 

"1 ,"2 

photons in one mode regardless of the number of photons 
in the other mode, Pa{n) is given by 



Pi(n) = ^P(ni,n2) 

^2H = E^K,ri2) 

jii 



(22) 



We can see easily that Pi{n) — P2{n) , a signature of the 
symmetric operation of the microlaser. In Fig. |3], we plot 
Pa (n) for Ub = 0.1, f = 50 and gnnt = 0.75, 0.8 and 1. 
Unlike the photon statistics of the single mode micro- 
laser, we notice in the probability distribution of a single 
mode of the two-mode microlaser. Pa (n) the existence of 
regions of exactly flat distribution. The explanation for 
these flat regions is plausible and goes as follows: 

The semi-classical rate equation at steady state of the 
three-level atom microlaser involves the total number of 
photons in the two modes (n+m) and reads: 



R sin^ [gTint V"- + 1 + + Ij = j[n + m] 



(23) 



Solutions of ( 23 ) specify steady state values for (n+m) . 



Since each value of the solutions for (n+m) can be formed 
by different combinations of n and m and all these com- 
binations are equally probable, it follows that the prob- 
ability distribution Pa (n) has a flat probability regions 
corresponding to all the different combinations of n and 
m. These flat regions are higher for lower n and steps 
down for larger values of n in a staircase pattern since 



FIG. 3: The photon number probability distribution 
Pa (n)for mode 1, calculated fornt — 0.1,^ — 50,gTint = 
0.75 (dashed), 0.8 (dotted), 1 (solid). 



for every possible solution for (n+m) the possible values 
of n and m start from zero while the high values of n and 
m are accessible only for large values of the solutions of 
(n+m) . 

It's notable to say that these flat probability sections 
are distinctive property of two- mode or higher microlaser. 
They are responsible for wider probability distribution 
than the poissonian distribution, i.e., super-poissonian 
distribution. In the more general case when the coupling 
between the atoms and the two modes are not identical, 
5i 7^ 52 , we will have the semi-classical rate equation 



, v/.gi (m+ 1) + .92 (n+ 1) = 7 [n + m] (24) 



It's evident that in this case no flat regions will appear 
in Pa (n) since we have a single n and m that satisfy 
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FIG. 4: The gain and loss parts of the semi-classical rate for 
grint = = 50,7 = 1 



the equation. This prediction will be confirmed shortly 
after applying the Quantum Trajectory Method to the 
microlaser problem. 



A. Second order correlation 

The steady state of the microlaser, as implictely men- 
tioned above, doesn't mean that the field inside the cav- 
ity of a certain microlaser setup has a fixed radiation 
intensity (i.e., a definite number of photons) since each 
atom will induce a change in the field and the field de- 
cay between two atoms induces further changes in the 
field or the number of photons. We are interested in this 
part to examine the correlation between these fluctua- 
tions in the cavity field, or more precisely, fluctuations in 
the number of photons inside the cavity represented by 
the second order correlation function g^-^^r). Measure- 
ment of the correlation function of the cavity field may 
provide a direct evidence of the quantized nature of light 
by detecting distinct correlation effects of the quantum 
field such as photon anti-bunching. The second order 
correlation function of the quantized field inside the cav- 
ity, (7^^'(t), is proportional to the probability of finding 
a pair of photons inside the cavity separated by a time 
T regardless of what happens to the photon number 
during this time interval and is defined by: 



. ^ _ (at(0)at(r)a(r)a(0)) 



(25) 



The correlation function of a micromaser cavity field 
can be calculated analytically as shown by Quang ^14^ 
by starting with a density matrix p(0), conditioned on 
the act of detecting and annihilating a photon from the 
cavity field. 



m 



apssO. 



where pss is the steady state solution of the master equa- 
tion. We evolve this conditional density matrix by the 
micromaser master equation and the correlation function 
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FIG. 5: The initial value of the second order correlation func- 
tion (j'^' (0) and the normalized mean photon number plotted 
versus the pump parameter 6 for a two-mode microlaser 



at time r will be proportional to the mean number of 
photons inside the cavity at time r as calculated by the 
evolved conditional density matrix pir) . 

(2j^^^ _ Tr [a^ap{T)\ 



V'V 

To calculate the evolution of the conditional density ma- 
trix p{t) , we solve the master equation by 4th order 
Runge Kutta Method. The quantum regression theorem 
has been used to derive this method. The other method 
to calculate g^'^^r) mimics what is done experimentally 
and is based on calculating the correlation between the 
times at which photons leak from the cavity. We will 
use the second method in the next section to calculate 
g'^^^r) numerically by the Quantum Trajectory Method. 

When we plot the initial value of the correlation func- 
tion of any of the two modes versus the pumping param- 
eter gTint , we notice as shown in Fig. |5|that g''^\0) is 
always higher than one except at the photon trapping 
states [IH] occurring at the severe dips in 5^^^(0) and in 
the normalized mean number of photons (n) . This means 
that photon number probability distribution (n) is 
super-poissonian for most of the range of gTi^t , a con- 
sequence of the flatness of Pa (n) explained earlier. This 
situation is not the same in the single mode microlaser 
due to the absence of this flatness. We depict in Fig. 
[6] the same graph for the single mode microlaser where 
the normalized mean number of photons and ^'^•'(O) are 
plotted versus the pumping parameter 6 deflned by: 



and 



gTint ' 



It's evident that the photon number distribution exihibits 
both sub-poissonian and super-poissonian statistics in 
the single mode microlaser. 

In the rest of this section we are going to concentrate 
on the correlation function at the trapping states deflned 
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FIG. 6: The initial value of the second order correlation func- 
tion (j*-^' (0) and the normalized mean photon number plotted 
versus the pump parameter 6 for a single-mode microlaser 
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FIG. 7: The second order correlation function (/'^''(r) corre- 
sponding to print = ^(b), ^(c)and^(d) calculated 
by QTM and compared with the analytical method 

by the condition gTi^t ^/n + 1 + rn + 1 = kir for a certain 
(n+m). At these states, the probabihty to find num- 
bers of photons larger than (n-|-m) will be identically 
zero since the probability for each atom to emit its pho- 
ton while interacting with the cavity field is proportional 
to sin^ [gTintViT' + 1 + w -I- l] . A trapping state is char- 
acterized by a sharp dip- called resonance- in the mean 
number of photons inside the microlaser cavity plotted 
versus the pumping parameter gTi^t ■ In Fig. [7] we show 
g^'^^r) calculated for the one, three and four photon trap- 
ping states and compared with the numerical calcula- 
tion by the Quantum Trajectory Method. Severe anti- 
bunching is noticed for the one-photon trapping state. 

We can understand the anti-bunching behavior of the 
one photon trapping state in terms of the necessary time 
needed between detecting a photon out from the cavity 
and re-pumping the cavity by an excited atom that de- 
posits another photon inside. This dead time between 
detecting a photon and re-pumping the cavity is respon- 
sible for the photons anti-bunching. This situation is very 
similar to the anti-bunching of the fiuorescence radiation 
emitted by a single atom where a dead time is unavoid- 
able between the emission of a photon and re-exciting 
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FIG. 8: The second order correlation function g^^\r) fitted 
with the analytical function /(r) 



the atom. We know from the theory of resonance fiuo- 
rescence of a single atom that the expression of the sec- 
ond order correlation g^'^^r) of the radiation scattered 
by a two-level atom driven by arbitrary field is given by 

mm 

(/'•^■'(t) = 1 — ^cos /IT + ^ sin /ir^ e~~'^ 

where F is the spectral linewidth of the atom, or alter- 
natively its spontaneous decay rate and fi is defined in 

terms of F and the Rabi frequency flfj hy n — \J ri|j — 

The analogy between the two-level atom and the one- 
photon trapping state is clear. The cavity plays the role 
of a two-level system where the two levels are either a 
photon is stored in the cavity or not. The difference be- 
tween the two systems is in the method of pumping and 
decaying from the higher level to the lower level. While 
the atom undergoes continuous Rabi oscillation by the 
pumping laser field, the cavity is pumped by a stream 
of atoms arriving randomly and separated by relatively 
large intervals. Probably this is the reason of not having 
the oscillatory behavior in the correlation function of the 
cavity field. We tried to fit g'^'^^ir) for this particular case 
with an analytical function and found that the function 
1 — e^'^'^ , where rj = R sin^ [gTintV^ +7 fits excellently 
with 5^^-'(t) as shown in Fig. [s] We found that the 
correlation function for the total number of photons (in 
the two modes) exhibits the same behavior and can be 
fitted with the same function. This behavior is not per- 
tinent to the two-mode microlaser, but appears also in 
the single-mode microlaser operating in its one photon 
trapping state characterized by sin^ [gTintV^ + l] = 0- 
We found that its correlation function fitted excellently 
with the analytical function 

/(r) = l-e-''\ (26) 

where rj = R siv? [gr-^^V^ + 7- We give a proof of this 
relation in the appendix 
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III. APPLYING THE QUANTUM 
TRAJECTORY METHOD TO SOLVE THE 
THREE-LEVEL ATOM MICROLASER 

A numerical method, basically a Monte Carlo simu- 
lation applied to quantum systems, to solve dissipative 
master equations was developed by three groups nearly 
at the same time in the early nineties [5J [HI UHl ITT] - 
In this method, called Quantum Trajectory Method 
(QTM), the observables of the system are obtained by 
averaging over many possible histories of the evolution 
of the system density matrix as a function of time. 
Each of these histories is called a trajectory, and its 
evolution is of a stochastic nature. Due to the statistical 
nature of quantum mechanics, taking the average over a 
large number of trajectories is equivalent to solving the 
master equation for this system. This is the essence of 
the quantum trajectory method. A certain trajectory 
can describe the stochastic evolution of the wavefunction 
or the density matrix of the open quantum system 
subjected to random quantum jumps representing its 
interaction with the reservoir. This method has a 
numerical advantage by reducing the computational 
power required to solve the master equation by the 
order of the dimension of the system, especially when 
the system has many degrees of freedom . Another 
advantage of the quantum trajectory method is the high 
level of control it allows on the parameters of the system. 
In our case we will use this ability to let each atom 
in any trajectory have a different velocity according to 
the velocity probability distribution, simulating what 
happens in reality. Including the variation of the atoms' 
speeds in the analytical solution of the master equation 
is a difficult task. Instead of evolving the density 



matrix in each trajectory, we will evolve wavefunctions 
representing the state of the cavity field. The simplest 
wavefunction one can use to represent the quantized 
electromagnetic field is the number state [TU]. For this 
reason and generalizing the quantum trajectory algo- 
rithm applied to the single mode micromaser developed 
by Pickles and Cresser in |10j , we evolve two number 
states |m) , |n) representing the deterministic number 
of photons in the two modes. In the most general case, 
where the two-mode micromaser cavity is maintained 
at a very low temperature, we can infer seven different 
events that may occur to the number state |m, n) with 
the following probabilities: 

1- An atom emits a photon in 
mode: |m, n) ^ |m-|-l,n) with 

(m -I- 1) sin^ [X{m, n)nr,t]/ ^■^{m, n) 

2- An atom emits a photon in 
mode: |m, n) — > |m, n -I- l)with 

l)gl sin^ [X{m,n)Tint]/ \'^{m,n) 

3- A photon from the first mode leaks out from the cavity 
\m,n) — + |m — with probability ^7i(".bi + l)m 

4- A photon from the second mode leaks out from 
the cavity |rn,n) — s- |m, n— 1), with probability 

^72(«fc2 + 1)"- 

5- A photon from the cavity walls is transferred to the 
cavity field of the first mode, \m,n) — > |m + l,?i) with 
probability ^7irtf,j (m -I- 1) 

6- A photon from the cavity walls is transferred to the 
cavity field of the second mode, |m, n) |m, n + I) with 
probabihty ^72^62 ('^ + 1) 

7- An atom passes through the cavity without emit- 
ting any photons, \m,n) — > \m,n), with probability 
^i?cos^ [A(m, n)rint] where 



the first 
probability: 

the second 
probability: 



R + '-fi [(nfcj -f- l)m -I- nt^ (m -t- 1)] + 72 [(f^h2 + 1)" + ^^62 ('^ + 1)] 

I 



The values of these probabilities can be derived by writ- 
ing the master equation in the Lindblad form 



and the waiting time between two successive jumps prob- 
ability distribution function is given by: 



m '- 

(27) 

where C™ represents the jump operator representing 
event (m). The probabihty of event (m) is calculated 
by: 



p(m) 



mr)\ClCMti)) 

Y.mi)\cl.c^\'<P{h)) 



(28) 



/ Y.{Ht)\clCr„\Ht))dt 



p{t) — e ' 
and for our system p{t) becomes , 



p{t) 



-At 



The choice of number states as the propagated wavefunc- 
tions has the advantage that the effective Hamiltonian 
controlling the evolution of the wavefunction between 
jumps keeps the number states unchanged [lU]. For more 
on the Quantum Trajectory Method, see [8l l9l ITT|. 

After generating many trajectories, we determine the 
diagonal elements of the density matrix by making a his- 
togram over the final states |m, n) of each trajectories. 
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FIG. 9: The photon number probability distribution 
Pa (n)for mode 1 calculated numerically and compared to the 
analytical solution for ii/7 = 200, gr — 0.3 




FIG. 10: The mean number of photons in mode 1 calculated 
numerically and compared to the analytical solution 



We show the photon number probability distribution cal- 
culated by the Quantum Trajectory Method for a micro- 
laser operating at R/^ — 200, gr — 0.3 in Fig. [9] where 
QTM shows an acceptable accuracy. 

The accuracy becomes better for lower values oi R/j . 
In Fig. 



10 we show the mean number of photons in one 
mode versus gTjnt obtained numerically and analytically. 
We used QTM to confirm that a microlaser operating in 
the non-symmetric mode, where the detailed balance con- 
dition doesn't apply, does reach steady state by checking 
that the avergae steady state density matrix calculated 
by QTM is the same when taking different lengthes of 
the trajectories. 

In fig. [TT]we see clearly that the flat regions disappear 
in the non-symmetric operation of the microlaser as pre- 
dicted earlier. It turns out that this flatness in Pa (n) 
is very sensitive to the difference between the coupling 
constants of the two modes as in flg. [12] 

To calculate the second order correlation function us- 
ing the Quantum Trajectory Method, we use a numeri- 
cal method similar to the one used experimentally |16j . 
In the experiment conducted by the group of M. Feld 
in MIT, the coherence function of the microlaser is ob- 
tained by calculating the correlation between the times 
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FIG. 11: The photon number probability distribution calcu- 
lated numerically for gi = g2 — 0.8 (solid) and gi = 0.8, 32 = 
0.5 (dashed) 




FIG. 12: The photon number probability distribution calcu- 
lated numerically for pi = (72 = 0.8 (solid) and gi = 0.8,52 = 
0.79 (dashed) 



when photons are emitted out from the cavity and de- 
tected by the photo-detector. Numerically, we have full 
details about each trajectory including the times at which 
photons leak from the cavity. So to calculate g'^'^^r), we 
compute the correlation between these times in each tra- 
jectory and average over all the trajectories. We have al- 
ready shown in Fig. [rjthe numerical calculation of g*^^^ (t) 
compared with the analytical calculation for the trapping 
states where g^^^r) exhibits anti-bunching behavior. We 
show in Fig. [13] the correlation function for two values of 
gTint w here the cavity field exhibits bunching behavior. 



In Fig. 



14 



we illustrate an interesting feature of g^'^^r) 
where the cavity field exhibits a transient anti-bunching 
behavior before it decays monotonically to one. 



IV. EFFECT OF ATOMIC VELOCITY 
DISTRIBUTION ON THE STATISTICAL 
PROPERTIES OF THE MICROLASER 

As we mentioned briefly in the introduction, atomic ve- 
locity selectors are not perfect and eventually the atoms 
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FIG. 13: The second order correlation function g^^\r) cal- 
culated for Ub = 0,^ = 10, QTint = 0.12 (a) and 0.5 (b) 
numerically and analytically. 
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FIG. 14: The second order correlation function g^^Hr) cal- 
culated for Ub — 0, ^ = 10, grint = 3 (a), 2 (b) and 1 (c) 
numerically and analytically. 



passing through the microlaser cavity will have some ve- 
locity distribution. In this section, we are going to il- 
lustrate the effect of this velocity spread of the atoms 
on the statistical properties of a single field mode of the 
microlaser cavity. It might be expected that the veloc- 
ity spread of the atoms will destroy the flat probability 
regions highlighted in the previous sections, but a quick 
look at Fig. [3] tells us that this is not correct. In fact, 
including a variety of interaction times is going to aver- 
age the flat regions in the probability distribution cor- 
responding to each value of ri„t and we end up with a 
persistent flat probability distribution whose width and 
height is a function of the relative atomic velocity spread 
as shown in Fig. |15| In this figure, we show that a rela- 
tive velocity spread of 20% maintains the flat regions in 
Pa (n), corresponding to gTint=0-Si and r/7 = 50. 

A vacuum trapping state is special trapping state 
where the cavity field is trapped at the vacuum state 
and occurs when the condition 5Ti„t\/l + 1 = n-T^ applies. 
This state, like other trapping states is characterized by 
a sharp dip in the plot of the mean number of photons. 
It is evident that the randomness in the interaction time 
Tint will destroy this condition and remove the resonances 
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FIG. 15: The photon number probability distribution calcu- 
lated numerically for the case of pumping the microlaser with 
a mono-velcoity atomic beam (solid) and an atomic beam 
having relative velocity spread of 20% (dotted) 
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FIG. 16: The mean number of photons at a trapping state 
for different relative velocity spreads. 



from the microlaser behavior. This is shown in Fig. 16 
where the mean number of photons in a vacuum trapping 
state is plotted for relative velocity distribution widths 
of .002%, .02%, .1%, 0.2%, 1% and 2%. It is notable that 
the trapping state is very sensitive to velocity distribu- 
tion width. 

We noticed the existence of sharp transitions in the 
mean number of photons plotted versus grint in Fig. |10[ 
These transitions occur between stationary solutions of 
the semi-slassical rate equation and we expect them to 
be induced by the randomness involved in the quantum 
system. We noticed, however, that including the spread 
in the atoms' velocities and hence increasing the ran- 
domness in the interaction times destroys these transi- 
tions starting from the transitions at large values of grint 
which are very sensitive to the velocity broadening. We 
show this behavior in Fig. [17] where the mean number 
of photons is plotted versus gTyat for a Guassian distri- 
bution of width 60%. This plots confirm that the system 
becomes more classical as more velocity fluctuations are 
introduced [17 . We can understand the immunity of 
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FIG. 17: The mean number of photons of the two-mode 
microlaser when relative velocity spreads of 60% is included 
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FIG. 18; The mean number of photons plotted versus the 
coupling constant g for different values of the interaction 
times to illustrate how a distribution of the interaction times 
(or equivalently atomic velocities) affect the microlaser phase 
transitions, tq is an arbitrary value of the interaction time 
equal to 1 s. The figure is for illustration and the values of g 
and To are not realistic. 

the first transition to the velocity distribution and the 
fragility of the peaks at higher values of grint by plotting 
the mean number of photons versus g for different values 
of the interaction times Ti^t- We see from the plots in 
Fig. [18] which correspond to a relative velocity spread of 
40%, why the first phase transition is not much affected 
while the higher transitions are easily destroyed. 

A. Effect of velocity spread on the correlation 
function 

We have seen in the theory of the three-level two- 
mode microlaser that the photon statistics of the cavity 
field exhibits bunching behavior for most of the range of 
the pumping parameter except at some of the trapping 
states, where the field is anti-bunched. We are going to 
show now the effect of the velocity spread on the second 
order coherence function in the two cases. We distinguish 
two regions from Fig. [Sj the first is the smooth region 
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FIG. 19: Second order correlation function for different 
atomic velocity spreads far from a trapping state — = 
0%, 20%, 50% and 100%. 

where g^^^ (0) < 2 and the second one is the sharp peaks of 
5^^-' (0) at the vacuum trapping states where gTint is multi- 
ple of tt/a/S. For the first case, we observed that the cor- 
relation is very immune to the atomic velocity spread and 
even a very broad velocity distribution doesn't reduce the 
field correlation substantially as shown in Fig. [T9| where 
the correlation due to mono-velocity atomic beam for a 
microlaser operating at grint = 0.6 and R = 10 is com- 
pared with microlaser pumped by an atomic beam having 
a spread of ^ = 0%, 20%, 50% and 100%. While in- 
creasing the velocity spread changes the average velocity 
for the atoms and hence drifts the operating point of the 
microlaser slightly, what we want to emphasize from Fig. 
[To] is that correlation is not affected much by the velocity 
spread. 

Finite correlation for practical atomic beams of rela- 
tive velocity spreads upto 20% has indeed been measured 
for the single-mode microlaser by Aljalal [16]. As for 
the regions near the vacuum trapping states, g^'^\0) is 
peaked because the number of photons inside the cavity 
is very small. The correlation function at these regions is 
strongly bunched which means that whenever a few num- 
ber of photons happen to exist inside the cavity, they will 
tend to leave the cavity together as a bunch of photons. 
At the extreme case of the vacuum trapping state, the 
correlation function diverges since the cavity has no pho- 
tons. We noticed that when operating the cavity near 
a vacuum trapping state where g'-^^^O) is very large, the 
correlation is very sensitive to the velocity spread of the 
atoms and collapses very fast until some residual corre- 
lation persists at a relative velocity spread of 2%. In 
Fig. |20] we show the second order correlation function 
of the cavity field g^^^r) for a microlaser operated close 
to this state {gunt = 1.9947r/A/3). In this figure, g^'^^r) 
of a mono-velocity beam is shown in addition to veloc- 
ity spreads of 0.02%, 0.03%, 0.04%, 0.06%, 0.08%, 
0.1%, 0.2% and 2%. We note that the correlation length 
has approximately the same value for first six values of 
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FIG. 20: The second order correlation function for different 
velocity spreads near a vacuum trapping state ^= 0.02%, 
0.03%, 0.04%, 0.06%, 0.08%, 0.1%, 0.2%, 2%. 



The reason for the fragihty of the correlation near 
a vacuum trapping state to the atomic velocity spread 
is that the smallest distribution in the interaction times 
destroys the vacuum trapping state as we saw in Fig. [16] 
and will introduce a small number of photons inside the 
cavity. These photons will cause the residual correlation 
mentioned above. Increasing the relative velocity spread 
beyond 2% doesn't affect this residual correlation consid- 
erably as in the first case (non-trapping states) where the 
correlation is very immune to the atomic velocity spread. 

On the other hand, when we investigated the effect 
of the atomic velocity spread on the correlation function 
for an operating point exhibiting anti-bunching, i.e., the 
one-photon trapping state, we found an interesting phe- 
nomenon. Adding more fluctuations in the atoms' ve- 
locity, produces a peak in the correlation function near 
T = and the field gradually becomes more correlated 
up to a relative velocity spread of ^ = 0.2% as shown 
in Fig. 21 In other words, the anti-bunching is con- 
verted into bunching. The explanation of this weird be- 
havior of noise-induced correlation goes as follows: when 
the velocity spread is slightly increased starting from the 
mono- velocity case, there will be a very small probability 
to find numbers of photons inside the cavity higher than 
one photon, i.e., a bunch of photons. Although the prob- 
ability is very small compared to the probability to find 
one photon or zero photon, its effect is overwhelming and 
eventually the correlation function is dominated by these 
rare bunches of photons that can exist inside the cavity. 

In a quantitative way, we can explain it from the initial 
value (?^^-*(0) which is given by [l^: 



(n^^ — {n) 



in) 



For the one photon trapping state, we have 
P(0),P(1) ^ and P{n) = for n > 1 and hence 
(n^) = (n) = P(l) This makes ff^^HO) = 0. For the 
slightest velocity distribution, the trapping state will be 
destroyed and P(n) will no longer be for n>l. This 

We can il- 
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makes ^^^•'(O) > as we see from Fig. 
lustrate this by a numerical example, for the case of a 



FIG. 21: The second order correlation function for a 
one-photon trapping state and relative velocity spreads 
Av/v = 0.01% (a), 0.03% (b), 0.05% (c),0.07% (d), 0.1% (e), 
0.15% (f) and 0.2% (g) 



relative velocity spread of .01%. We find in this case 
from the QTM simulation that P{n) > for n < 7 and 
(n) = 0.3918 and (n^) = 0.4264. These values yield 
an initial value of the correlation function g'^'^^Q) to be 
0.2254. It turns out that the wider the velocity distribu- 
tion, the higher the value of ^''^-'(O) is till Av/v reaches 
0.2%. By increasing the relative velocity spread above 
0.2%, the correlation is lost gradually due to the huge 
randomness in the interaction times between the atoms 
and the cavity, until a residual correlation persists start- 
ing from a relative velocity spread of 20%. This is clear 
in Fig. |22] where we notice that correlation and hence 
the bunching of the cavity photons decreases gradually 
for relative velocity spreads of 0.2%, 0.6%, 1%, 2% and 
20%. Next, we show the effect of velocity distribution on 
the correlation function of the total number of photons 
in the cavity (i.e., in the two modes combined). With a 
velocity spread of 0.4%, and gTiat = -% and we get 
the plots shown in Fig. |23] where as in Fig. |2T} we no- 
tice that the anti-bunching behavior in .g'-^-'(T) has been 
converted into bunching for the same reason mentioned 
above. 



V. CONCLUSION 

We have applied the Quantum Trajectory Method to 
the two mode microlaser operating on A-type three level 
atoms. We verified that the two-mode microlaser does 
reach steady state when the coupling between the atom 
and the two modes is not symmetric, the case where the 
detailed balance condition doesn't apply. As for the sym- 
metric operation of the microlaser, we explained the flat 
probability regions in the photon number probability dis- 
tribution of any of the two modes and emphasized the 
fact that the existence of two modes equally coupled to 
the atoms gives rise to these flat regions. The super- 
poissonian distribution of the photon statistics of any of 
the two modes is a direct consequence of these flat re- 
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FIG. 22: The second order correlation function for a 
one-photon trapping state and relative velocity spreads 
Av/v = 0.2%(a), 0.6%(b), l%(c), 2%(d)and20%(e)The micro- 
laser is operated at print ~ 1.81379, R — 10, n^= 




FIG. 23: The second order correlation function for the total 
number of photons for the one-photon trapping states grint = 
tt/VS (thick) and 27r/\/3 (thin) and relative velocity spreads 
Av/v = 0% (a) and 0.04% (b) 



gions and causes the photon correlation to be bunched 
for most of the range of the operating parameters of the 
microlaser, a distinct property that distinguishes the two- 
mode microlaser from the single mode microlaser. We in- 
troduced the analogy between the correlation of the field 
of the one-photon trapping state and the resonance flu- 
orescence radiation scattered from a two-level atom. We 
fitted the correlation function of this state to a very sim- 
ple formula and give an analytical derivation of this for- 
mula in the appendix. Experimental measurement of this 
anti-bunching is a direct verification of the quantized na- 
ture of the field. We investigated the effect of the atomic 
velocity spread on the photon statistics by the Quantum 
Trajectory Method. We found that the trapping states 
are very sensitive to the velocity spread of the atoms. 
The sharp transitions in the microlaser output are also 
destroyed by the velocity distribution of the atoms start- 
ing from the transitions at higher values of the pumping 
parameter. The second order correlation is in general 
immune to the velocity spread of the atoms. For the 
special case when the microlaser is operating near a vac- 



uum trapping state, the intensity correlation of the cav- 
ity field is easily destroyed by a relative velocity spread 
as low as 0.2 %. Increasing the velocity spread further 
doesn't affect the low residual correlation that survives 
the velocity spread of the atoms. When the cavity field of 
the mono-velocity atomic beam microlaser exhibits anti- 
correlation, adding randomness in the atomic velocities 
counter-intuitively converts the anti-correlation into cor- 
relation. This noise-induced correlation stops at velocity 
spreads of about 0.2 %. Increasing the velocity spread 
further destroys the correlation very slowly until a resid- 
ual correlation persists for relative velocity spreads of 
20% and higher. 

We propose as a further investigation developing a 
quantum trajectory method from the Fokker-Planck 
equation of the microlaser |17j and evolving coherent 
states of radiation instead of number states, as the coher- 
ent states are the closest states to the classical radiation 
field. 

The quantized motion of the atom inside the cavity 
could possibly be included in the numerical simulation by 
quantizing the longitudinal dimension of the cavity and 
involving the position of the atom in the wavefunction 
we propagate by QTM during the atom field interaction. 
The quantized motion of the atom inside the cavity was 
shown to lead to some distinct features of the two-mode 
micromaser jl8) . 

Lastly, we would like to propose a novel numerical 
method to solve the master equation of the microlaser 
by using a Genetic Algorithm. The basic idea behind 
using Genetic Algorithms to solve complex problems is 
to generate an initial random ensemble of solutions and 
then apply some criterion or fitness condition to select 
some of them that are closer to the proper solution of the 
problem. These good candidates are merged together, 
swapped, paired, or mutated to form another generation 
of solutions upon which we apply the same fitness cri- 
terion to select the good candidates and repeat this cy- 
cle till we obtain the best fitting solution to the prob- 
lem. When applying this algorithm to the microlaser, 
the elements of the ensemble of solutions will be the den- 
sity matrices representing the steady state field statistics. 
The criterion would be minimizing the change incurred 
to the matrix when an atom passes through the cavity 
and the field is let to decay. The ideal solution should 
have zero change from the definition of the micromaser 
steady state. We select the good candidates of any gen- 
eration of solutions by picking the matrices that yield the 
minimum of this change. Merging a set of matrices to- 
gether to bring up another generation of solutions should 
be an easy numerical task. 



APPENDIX 



We are going to prove that for the one photon t rapp ing 



state, '{t) is exactly equal to /(r) given in (26 1 for 



the case of the two level single mode microlaser since it's 
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much easier than the two- mode case. Starting from the 
definition of g^^^r) for a quantized field: 



(at(0)at(r)a(r)a(0)) 



(A.l) 



where r = represents the steady state. For the one pho- 
ton trapping state only the states |0) , |1) are accessible 
to the cavity field. We can then write 

E {n\aHO)aHT)a{T)a{0)\n) 

(2)/ \ n=0,l 

J2 n. {n — 1\ a^{T)a{T) \n — 1) 

n=0,l 



|at(r)a(r) 



\nir)\ 



(A.2) 



where n(r) is the number operator a^{T)a{T). Under 
the one-photon trapping state condition, we have only 
two possibilities: to find or 1 photon inside the cavity 
and hence 



i{t) = ^n.pn{t) ^Pi{t) 



(A.3) 



Since pi{t) -t-po(^) = 1, we can write (A. 5 1 as: 



i{t) = (1 - n{t))Ri 



i?sin^ gTintV2 ~ (^Rsin^ gTintV2 + n{t) 



where 



\ntV2 



gnnt 



A - Bn{t) 



(A.7) 



int\/2 



By integrating (A.7 1 from ^ t we get 



dn 



A-B.n 

n{0) 

By solving for n^r), we obtain: 
'A 



dt 



( \ ^ 
nKT) = 

Since we know that 



B 



- "(0) 



i(0) = ^(n|n(0) 



+ 7 



(A., 



(A.9) 



(A.IO) 



(A.ll) 



since Pn{i) = for n ^ 0, 1 . It can be shown from the 
single mode microlaser master equation (see for example 
[13]) that evolution of pi{t) is governed by: 



gnnt 



(A.4) 



From ( A.4[ A.3 1, we deduce that the evolution of the 
mean number of photons inside the cavity is governed by 



i(t) = ^ {n\n{T) 



(A. 12) 



n(0) ^ ^ {n\ n(0) \n) \n) (n| (A.13) 

n 

and by assuming that n(r) is diagonal in the form 



n{t) — po{t)Rs'm^ 



.gTint 



(A.5) 



This equation is intuitive for the one-photon trapping 
state and could have been written directly without re- 
ferring to master equation. The second term on the left 
hand side represents the rate of photon loss from the cav- 
ity while the first term represents the number of photons 
injected inside the cavity per unit time. At steady state, 
pi(0) = and hence 

= {1 -pi(0)}i?sin2 \gr.^^^] -7pi(0) 



which leads to 



(n) = (n2)=pi(0) = 



flsin^ [ffTi„t\/2] 
7 + i?sin^ [gri„tV2 



(A.6) 



n(r) = ^ {n\ h^r) \n) \n) (n| (A.14) 

n 

,we can infer from ( |A.10|A.ll|A.12|A.13|A.T4l ) that 



^ ' B 



-Bt 



(A.15) 



and hence the numerator of (A.2 1 equals 
'A 



\n{r)\0) = ^- 



B 



- 



-B 



A 



(A.16) 



From ( |A.16[ ), ( |A.2[ ), ( |A.6[ ) and dO] ) we can write 

g(2)(r)=(l-e-«0 (A.17) 
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